from kuai.molio import read_mol, write_mol
from kuai.ffio import read_ff, write_ff

import os
import os.path

DFF_FREQ_DFI = """$SYSTEM
   ROOT = %(DFF_ROOT)s
$END
$MODEL
   %(TEMP_FILE_NAME)s.msd
$END
$FORCEFIELD
   TYPE = AMBER
   INPUT = %(TEMP_FILE_NAME)s.ppf
   USE_CUTOFF = NO
   EQUI = INPUT
$END
$COMPUTE
  CALC = FULL
$END
$PRINT
   COORDINATE = NO
   FREQUENCE = YES
$END
"""

def hessian(mol, ff, vars = {}):
    if "DFF_BIN" not in vars:
        vars["DFF_BIN"] = r"F:\Madlee\Job\Direct Force Field\7\bin32w"
    if "DFF_ROOT" not in vars:
        vars["DFF_ROOT"] = os.path.dirname(vars["DFF_BIN"])
    if "TEMP_FILE_NAME" not in vars:
        vars["TEMP_FILE_NAME"] = ".dff"

    dfi = open(vars["TEMP_FILE_NAME"] + ".dfi", "w")
    dfi.write(DFF_FREQ_DFI % vars)
    dfi.close()
    write_mol(vars["TEMP_FILE_NAME"] + ".msd", mol)
    # write_ff(vars["TEMP_FILE_NAME"] + ".ppf", ff)

    command = '"' + os.path.join(vars["DFF_BIN"], "dffmm") + '" '
    command +=  " " + vars["TEMP_FILE_NAME"]
    pin, pout, perr = os.popen3(command)
    pin.close()
    pout.read()
    error = perr.read().strip()
    if error != "":
        raise RuntimeError(error)
    else:
        output = open(vars["TEMP_FILE_NAME"] + ".dfo")
        result = []
        for i in output:
            if i.startswith("The second derivatives"):
                for j in output:
                    if j.strip() == '':
                        break
                    else:
                        result += [float(k.strip()) for k in j.split()]
                assert len(result) == len(mol.atoms) * 3 * (len(mol.atoms) * 3 + 1) / 2
                return result

        if len(result) == 0:
            raise RuntimeError('There are no hessian data in the result')
        else:
            return result
        return read_mol(vars["TEMP_FILE_NAME"] + ".msd")

